UAS Presentation for Rio Chama CFLRP
George Woolsey
21 November, 2024
New Mexico Example
Let’s go through the tree clump grouping process for a single stand flown on the Santa Fe National Forest in November 2024.
Setup
# library
library(tidyverse) # the tidyverse
library(viridis) # viridis colors
library(scales) # work with number and plot scales
library(latex2exp) # math formulas text
# library(paletteer) # working with color palettes
# library(palettetown) # color palettes
# palettetown::pokedex(cb = 4)
library(harrypotter) # color palettes
# spatial analysis
library(terra) # raster
library(sf) # simple features
# visualization
library(kableExtra)
library(mapview) # interactive html maps
library(patchwork) # ! better align plots in grid
library(ggnewscale) # add another scale to ggplot
library(ggtext) # customize plot text colors
library(ggpubr) # alternative plot fn's (table)
# forest patches
# install.packages("pak")
# pak::pkg_install("mhahsler/dbscan", upgrade = T)
library(dbscan)now, define a function to convert metric to imperial units
calc_imperial_units_fn <- function(df) {
df %>%
# convert to imperial units
dplyr::mutate(
dplyr::across(
.cols = tidyselect::ends_with("_cm")
, ~ .x * 0.394
, .names = "{.col}_in"
)
, dplyr::across(
.cols = tidyselect::ends_with("_m")
, ~ .x * 3.28
, .names = "{.col}_ft"
)
, dplyr::across(
.cols = tidyselect::ends_with("_m2_per_ha")
, ~ .x * 4.359
, .names = "{.col}_ftac"
)
, dplyr::across(
.cols = tidyselect::ends_with("_per_ha") & !tidyselect::ends_with("_m2_per_ha")
, ~ .x * 0.405
, .names = "{.col}_ac"
)
, dplyr::across(
.cols = tidyselect::ends_with("_area_ha")
, ~ .x * 2.471
, .names = "{.col}_ac"
)
, dplyr::across(
.cols = tidyselect::ends_with("_m2")
, ~ .x * 10.764
, .names = "{.col}_ft2"
)
) %>%
dplyr::rename_with(
.fn = function(x){dplyr::case_when(
stringr::str_ends(x,"_cm_in") ~ stringr::str_replace(x,"_cm_in","_in")
, stringr::str_ends(x,"_m_ft") ~ stringr::str_replace(x,"_m_ft","_ft")
, stringr::str_ends(x,"_m2_per_ha_ftac") ~ stringr::str_replace(x,"_m2_per_ha_ftac","_ft2_per_ac")
, stringr::str_ends(x,"_per_ha_ac") ~ stringr::str_replace(x,"_per_ha_ac","_per_ac")
, stringr::str_ends(x,"_area_ha_ac") ~ stringr::str_replace(x,"_area_ha_ac","_area_ac")
, stringr::str_ends(x,"_m2_ft2") ~ stringr::str_replace(x,"_m2_ft2","_ft2")
, TRUE ~ x
)}
)
}Load Data
load data
# where is the processed data from cloud2trees?
input_dir <- "../data/"
# tree top points
treetops_sf_with_dbh <- sf::st_read(paste0(input_dir, "/final_detected_tree_tops.gpkg"))## Reading layer `final_detected_tree_tops' from data source
## `C:\Data\usfs\rio_chama_cflrp_202411\data\final_detected_tree_tops.gpkg'
## using driver `GPKG'
## Simple feature collection with 80693 features and 17 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 341046.3 ymin: 4006169 xmax: 342265.5 ymax: 4007337
## Projected CRS: NAD83 / UTM zone 13N# what?
treetops_sf_with_dbh %>% dplyr::glimpse()## Rows: 80,693
## Columns: 18
## $ treeID <chr> "1_341396.7_4007336.9", "2_341418.7_4007336.…
## $ tree_height_m <dbl> 14.132, 8.233, 10.531, 2.682, 5.747, 3.252, …
## $ crown_area_m2 <dbl> 1.16, 1.68, 0.56, 0.16, 0.36, 0.16, 0.16, 0.…
## $ fia_est_dbh_cm <dbl> 32.854644, 15.114067, 20.838229, 4.806052, 9…
## $ fia_est_dbh_cm_lower <dbl> 13.800244, 6.434515, 8.645597, 1.980668, 4.0…
## $ fia_est_dbh_cm_upper <dbl> 58.496900, 26.748053, 37.532606, 8.588588, 1…
## $ dbh_cm <dbl> 32.854644, 15.114067, 20.838229, 4.806052, 9…
## $ is_training_data <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ dbh_m <dbl> 0.32854644, 0.15114067, 0.20838229, 0.048060…
## $ radius_m <dbl> 0.16427322, 0.07557033, 0.10419115, 0.024030…
## $ basal_area_m2 <dbl> 0.084778049, 0.017941244, 0.034104486, 0.001…
## $ basal_area_ft2 <dbl> 0.91255092, 0.19311955, 0.36710069, 0.019527…
## $ ptcld_extracted_dbh_cm <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ptcld_predicted_dbh_cm <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_trees_per_ha <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_relative_tree_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_dist_to_nearest_m <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ geom <POINT [m]> POINT (341396.7 4007337), POINT (34141…# where?
treetops_sf_with_dbh %>%
dplyr::slice_sample(prop = 0.01) %>%
mapview::mapview()create a fake stand boundary
# let's just zoom in on the center ha and call this our "stand"
stand_sf <- treetops_sf_with_dbh %>%
sf::st_bbox() %>%
sf::st_as_sfc() %>%
sf::st_centroid() %>%
sf::st_buffer(dist = sqrt(20000)/2, endCapStyle = "SQUARE")
# area
sf::st_area(stand_sf)## 20000 [m^2]check this stand on the map
zoom out if the map looks blank
mapview::mapviewOptions(
homebutton = FALSE
, basemaps = c("Esri.WorldImagery","OpenStreetMap")
)
mapview::mapview(stand_sf, alpha.regions = 0, color = "red", lwd = 3)zoom out if the map looks blank
let’s make a base plot of the stand
plt_stand <- ggplot() +
geom_sf(data = stand_sf, color = "black", fill = NA) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(
x = ""
, y = ""
) +
theme_void() +
theme(
legend.position = "top" # c(0.5,1)
, legend.direction = "horizontal"
, legend.margin = margin(0,0,0,0)
, legend.text = element_text(size = 8)
, legend.title = element_text(size = 8)
, legend.key = element_rect(fill = NA, color = NA)
# , plot.title = ggtext::element_markdown(size = 10, hjust = 0.5)
, plot.title = element_text(size = 10, hjust = 0.5, face = "bold")
, plot.subtitle = element_text(size = 8, hjust = 0.5, face = "italic")
)let’s define our DBH limit for classifying trees as “overstory” and we’ll define our distance in meters to separate tree clump groupings
tree_clump_dist_m <- 6
ostory_dbh_cm <- 6*2.54 # = inches*2.54plot the stand with some trees
plt_stand +
geom_sf(
data = treetops_sf_with_dbh %>%
sf::st_intersection(stand_sf %>% sf::st_buffer(10)) %>%
# overstory
dplyr::mutate(overstory = dbh_cm>=ostory_dbh_cm)
, mapping = aes(color = overstory)
)
ggsave("../data/01_tree_tops.jpg", width = 7, height = 7)Define Clumps
create a general function to apply dbscan::dbscan to
point data of class sf
With respect to clump size groupings, Churchill et al. (2016) note that:
Proportions for clump sizes should be lumped into four or five bins for operational simplicity. We use 4 or 5 bins (Fig 5): individual trees, small clumps (2-4 trees), medium clumps (5-9 trees), and large clumps (10-20+ trees). Note that when instructed to leave a large clump (e.g. 10-20 trees), marking crews often have difficulty leaving the upper end of the size range (e.g. an 18, 19, or 20 tree clump). Thus adding a fifth bin for “super clumps” may be necessary (e.g. 15-20 trees or 20-25+ trees), especially if the upper size range of clumps is desired. (p. 12-13)
In their prescription worksheet Excel file Churchill et al. (2016) have the following group sizes: 1, 2-4 (Mean clump size: 3), 5-9 (7), 10-15 (12), 16-25 (20). To be able to cut clumps down to achieve the desired largest “16-25” group, we’ll add a “super” group with >25 trees. This same principle can be applied if you desire a clump size group larger than this (i.e. 25-35), whereby a group above this should be added so that the resulting prescription cuts the “super” group to ensure the proper sizing of the desired largest group.
# function to clump data that are sf points
st_clump_points <- function(
x # point data of class `sf`
, clump_dist_m = 6 # size (radius) of the epsilon neighborhood = maximum distance between points to add to clump
, clump_breaks = c(0,1,4,9,15,25,Inf) # where to break the clump size groups
, clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees") # what to name the clump size groups
) {
# get points as x,y
point_clusters = x %>%
dplyr::mutate(
X = sf::st_coordinates(.)[,1] %>% as.numeric()
, Y = sf::st_coordinates(.)[,2] %>% as.numeric()
)
#############################################################################
##### Identify clusters in each stem map plot #####
#############################################################################
### Place trees into clusters using an inter-tree distance of 6 m
my_dbscan_temp = point_clusters %>%
sf::st_drop_geometry() %>%
dplyr::select(X,Y) %>%
dbscan::dbscan(eps = clump_dist_m, minPts = 2)
# my_dbscan_temp %>% str()
### append cluster ID to tree points
point_clusters$dbscan_cluster = my_dbscan_temp$cluster
# point_clusters$cluster %>% summary()
# point_clusters %>% sf::st_drop_geometry() %>% dplyr::count(cluster) %>% dplyr::arrange(desc(n)) %>% dplyr::slice_head(n=11)
### cluster metrics
point_clusters = point_clusters %>%
dplyr::group_by(dbscan_cluster) %>%
dplyr::mutate(
# unique dbscan_cluster for individuals
clump_id = dplyr::case_when(
dbscan_cluster == 0 ~ max(my_dbscan_temp$cluster)+dplyr::row_number()
, T ~ dbscan_cluster
) %>%
factor()
) %>%
dplyr::group_by(clump_id) %>%
dplyr::mutate(
dbscan_cluster = factor(dbscan_cluster)
, clump_n_trees = dplyr::n()
, clump_n_trees_grp = cut(
clump_n_trees
, breaks = clump_breaks
, labels = clump_labels
) %>%
factor(
ordered = T
, levels = clump_labels
)
) %>%
dplyr::ungroup() %>%
dplyr::mutate(tree_clump_dist_m = clump_dist_m)
# return
return(point_clusters)
}function specific to this project to filter the tree points attached
to harvest bounds and apply the st_clump_points
function
# create function to pass a unit id and return list of trees with clump groupings
get_tree_clumps = function(
tree_clump_dist_m=6
, ostory_ht_m = as.numeric(NA)
, ostory_dbh_cm = as.numeric(NA)
){
# check ostory definition
if(is.na(as.numeric(ostory_dbh_cm)) & is.na(as.numeric(ostory_ht_m))){
warning("`ostory_dbh_cm` and `ostory_ht_m` are not set...using `ostory_dbh_cm` = 12.7")
ostory_dbh_cm = 5*2.54
# filter data
ttops_temp = treetops_sf_with_dbh %>% sf::st_intersection(stand_sf) %>%
dplyr::filter(
dbh_cm>=as.numeric(ostory_dbh_cm)
)
}else if(!is.na(as.numeric(ostory_dbh_cm))){
# filter data
ttops_temp = treetops_sf_with_dbh %>% sf::st_intersection(stand_sf) %>%
dplyr::filter(
dbh_cm>=as.numeric(ostory_dbh_cm)
)
}else{
# filter data
ttops_temp = treetops_sf_with_dbh %>% sf::st_intersection(stand_sf) %>%
dplyr::filter(
tree_height_m>=as.numeric(ostory_ht_m)
)
}
### Place trees into clusters using an inter-tree distance of 6 m
ttops_temp = st_clump_points(
x = ttops_temp
, clump_dist_m = tree_clump_dist_m
)
# add distance to nearest within clump
ttops_temp =
ttops_temp %>%
dplyr::group_by(clump_id) %>%
tidyr::nest() %>%
dplyr::mutate(
distance_clump_nn_m = purrr::map(data, function(x){
# get index of nearest neighbor
i = sf::st_nearest_feature(x)
# get dist
d = sf::st_distance(x, x[i,], by_element=TRUE) %>% as.numeric()
return(d)
})
) %>%
tidyr::unnest(cols = c(data, distance_clump_nn_m)) %>%
sf::st_set_geometry("geom") %>% # set it cuz it got lost
dplyr::ungroup() %>%
dplyr::mutate(
tree_clump_dist_m = tree_clump_dist_m
, clump_id_duplicate = clump_id # can use this even after nesting data by clump_id
# , ostory_ht_m = ifelse(is.na(ostory_ht_m), as.numeric(NA), as.numeric(ostory_ht_m))
# , ostory_dbh_cm = ifelse(is.na(ostory_dbh_cm), as.numeric(NA), as.numeric(ostory_dbh_cm))
)
# return
return(ttops_temp)
}call that function
# call it
ttops_temp = get_tree_clumps(
tree_clump_dist_m = 6
, ostory_dbh_cm = ostory_dbh_cm # 6 in dbh
)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometrieswhat?
ttops_temp %>% dplyr::glimpse()## Rows: 769
## Columns: 27
## $ clump_id <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ treeID <chr> "49736_341606.7_4006823.3", "49892_341603.5_…
## $ tree_height_m <dbl> 20.767, 9.622, 15.808, 12.870, 11.495, 10.35…
## $ crown_area_m2 <dbl> 32.72, 1.52, 6.28, 4.28, 3.52, 1.68, 8.48, 2…
## $ fia_est_dbh_cm <dbl> 61.52282, 18.49889, 38.74657, 28.38951, 23.8…
## $ fia_est_dbh_cm_lower <dbl> 25.860036, 7.815217, 16.136208, 12.057785, 9…
## $ fia_est_dbh_cm_upper <dbl> 109.37700, 32.93805, 68.86084, 51.30199, 42.…
## $ dbh_cm <dbl> 61.52282, 18.49889, 38.74657, 28.38951, 23.8…
## $ is_training_data <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ dbh_m <dbl> 0.6152282, 0.1849889, 0.3874657, 0.2838951, …
## $ radius_m <dbl> 0.30761410, 0.09249444, 0.19373285, 0.141947…
## $ basal_area_m2 <dbl> 0.29727772, 0.02687702, 0.11791157, 0.063300…
## $ basal_area_ft2 <dbl> 3.1998974, 0.2893042, 1.2692001, 0.6813641, …
## $ ptcld_extracted_dbh_cm <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ptcld_predicted_dbh_cm <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_trees_per_ha <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_relative_tree_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_dist_to_nearest_m <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ geom <POINT [m]> POINT (341606.7 4006823), POINT (34160…
## $ X <dbl> 341606.7, 341603.5, 341609.9, 341616.9, 3416…
## $ Y <dbl> 4006823, 4006822, 4006822, 4006822, 4006821,…
## $ dbscan_cluster <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ clump_n_trees <int> 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, …
## $ clump_n_trees_grp <ord> >25 trees, >25 trees, >25 trees, >25 trees, …
## $ tree_clump_dist_m <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,…
## $ distance_clump_nn_m <dbl> 3.352611, 2.154066, 3.352611, 2.630589, 2.72…
## $ clump_id_duplicate <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…Clump Polygons and Metrics
Churchill et al. (2016) provide instructions for implementing the clump identification (Plotkin et al. 2002) in ArcGIS:
Use the Buffer tool (in the Proximity toolset within the Analysis toolbox) to create a buffer of distance d/2, one half the inter-tree distance, around each point. This quantity d/2 is meant to approximate the crown radius of a “typical” overstory tree. Set the Dissolve Type option to ALL, which dissolves overlapping buffers, creating a reduced set of spatially non-overlapping polygons stored as a multipart polygon feature…Sanchez Meador et al. (2011) provide some useful examples of how clump attributes can be summarized…The method described here can be modified to use measured or modeled crown radii for each tree in place of d/2 (p.36)
# create function to pass a return from st_clump_points/get_tree_clumps and create clump polygons with summary stats
get_clump_summary = function(dta){
# get tree_clump_dist_m
tree_clump_dist_m = min(dta$tree_clump_dist_m, na.rm = T)
# create clump polys and summary
clump_polys_temp =
dta %>%
dplyr::ungroup() %>%
sf::st_set_geometry("geometry") %>%
sf::st_buffer(tree_clump_dist_m/2) %>%
dplyr::group_by(clump_id, dbscan_cluster, clump_n_trees_grp) %>%
dplyr::summarise(
# union buffered tree points
geometry = sf::st_union(geometry)
# summary metrics
, n_trees = dplyr::n_distinct(treeID)
, mean_dbh_cm = mean(dbh_cm, na.rm = T)
, mean_tree_height_m = mean(tree_height_m, na.rm = T)
, loreys_height_m = sum(basal_area_m2*tree_height_m, na.rm = T) / sum(basal_area_m2, na.rm = T)
, basal_area_m2 = sum(basal_area_m2, na.rm = T)
, sum_dbh_cm_sq = sum(dbh_cm^2, na.rm = T)
, .groups = "drop_last"
) %>%
dplyr::ungroup() %>%
sf::st_make_valid() %>%
dplyr::mutate(
clump_area_ha = sf::st_area(geometry) %>% as.numeric() %>% `/`(10000)
, trees_per_ha = (n_trees/clump_area_ha)
, basal_area_m2_per_ha = (basal_area_m2/clump_area_ha)
, pct_stand_basal_area = basal_area_m2/sum(basal_area_m2)
, pct_stand_n_trees = n_trees/sum(n_trees)
, qmd_cm = sqrt(sum_dbh_cm_sq/n_trees)
) %>%
dplyr::select(-c(sum_dbh_cm_sq)) %>%
# convert to imperial units
calc_imperial_units_fn() %>%
dplyr::mutate(tree_clump_dist_m = tree_clump_dist_m)
# calculate distance between clumps
clump_polys_temp = clump_polys_temp %>%
dplyr::mutate(
nearest = sf::st_nearest_feature(clump_polys_temp)
, distance_nearest_clump_m = sf::st_distance(
clump_polys_temp
, clump_polys_temp[nearest,]
, by_element=TRUE
) %>%
as.numeric()
) %>%
dplyr::select(-c(nearest))
# return
return(clump_polys_temp)
}get the clump summary data
# get it
# get_clump_summary(
# dta = get_tree_clumps(tree_clump_dist_m = tree_clump_dist_m)
# )
clump_polys_temp = get_clump_summary(ttops_temp)
# what?
clump_polys_temp %>% dplyr::glimpse()## Rows: 39
## Columns: 25
## $ clump_id <fct> 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 3…
## $ dbscan_cluster <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, …
## $ clump_n_trees_grp <ord> Individual, Individual, Individual, Individua…
## $ geometry <GEOMETRY [m]> POLYGON ((341593.9 4006821,..., POLY…
## $ n_trees <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 31, 83, 1…
## $ mean_dbh_cm <dbl> 96.40704, 35.15035, 94.58870, 97.47311, 66.14…
## $ mean_tree_height_m <dbl> 26.72900, 14.89200, 26.50900, 26.99600, 21.65…
## $ loreys_height_m <dbl> 26.729000, 14.892000, 26.509001, 26.996000, 2…
## $ basal_area_m2 <dbl> 0.72997403, 0.09703962, 0.70269748, 0.7462073…
## $ clump_area_ha <dbl> 0.002826142, 0.002826142, 0.002826142, 0.0028…
## $ trees_per_ha <dbl> 353.8393, 353.8393, 353.8393, 353.8393, 353.8…
## $ basal_area_m2_per_ha <dbl> 258.29350, 34.33643, 248.64199, 264.03750, 12…
## $ pct_stand_basal_area <dbl> 0.0071688845, 0.0009530008, 0.0069010087, 0.0…
## $ pct_stand_n_trees <dbl> 0.00130039, 0.00130039, 0.00130039, 0.0013003…
## $ qmd_cm <dbl> 96.40704, 35.15035, 94.58870, 97.47311, 66.14…
## $ mean_dbh_in <dbl> 37.984375, 13.849236, 37.267948, 38.404406, 2…
## $ qmd_in <dbl> 37.984375, 13.849236, 37.267948, 38.404406, 2…
## $ mean_tree_height_ft <dbl> 87.67112, 48.84576, 86.94952, 88.54688, 71.03…
## $ loreys_height_ft <dbl> 87.67112, 48.84576, 86.94952, 88.54688, 71.03…
## $ basal_area_ft2_per_ac <dbl> 1125.90139, 149.67250, 1083.83043, 1150.93946…
## $ trees_per_ac <dbl> 143.3049, 143.3049, 143.3049, 143.3049, 143.3…
## $ clump_area_ac <dbl> 0.006983396, 0.006983396, 0.006983396, 0.0069…
## $ basal_area_ft2 <dbl> 7.8574404, 1.0445344, 7.5638357, 8.0321761, 3…
## $ tree_clump_dist_m <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, …
## $ distance_nearest_clump_m <dbl> 2.12332457, 0.04846820, 0.01763506, 0.8018667…# do these numbers match
identical(
# clump
nrow(clump_polys_temp)
# clumps in tree list data
, ttops_temp %>% dplyr::distinct(clump_id) %>% nrow()
)## [1] TRUEplot the clump group size
plt_stand +
geom_sf(
data = clump_polys_temp
, mapping = aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
labs(fill = "") +
guides(size = "none", fill = guide_legend(override.aes = list(size = 5)))
ggsave("../data/02_tree_clumps.jpg", width = 7, height = 7)Clump Spacing
See Churchill et al. (2016) Figure 4 (p.10) and Matonis and Binkley (2018) who “calculated coverage of mosaic-meadows (percentage of stand > 6 m from overstory trees)” (p. 124)
Since we already buffered the tree points to approximate the crown radius, we’ll continue to use our \(d/2\) where \(d\) is maximum distance between trees for determining tree clumps and is meant to approximate the crown radius of a “typical” overstory tree
# create function to pass a return from get_clump_summary() and get a distance raster
get_clump_dist_rast = function(dta){
# get tree_clump_dist_m
tree_clump_dist_m = min(dta$tree_clump_dist_m, na.rm = T)
# rasterize the clump polygons and then calculate distance between clumps as raster
dist_rast =
terra::rasterize(
x = dta %>% terra::vect()
, y = dta %>%
terra::vect() %>%
terra::rast(res = 0.2)
) %>%
terra::distance() %>%
# crop it to stand extent
terra::crop(
stand_sf %>%
terra::vect()
) %>%
terra::mask(
stand_sf %>%
terra::vect()
)
######### part 2
# now create openings vector data
openings_vect =
dist_rast %>%
terra::classify(rcl = c(tree_clump_dist_m/2,Inf), others = NA, include.lowest = T) %>%
terra::as.polygons(na.rm = T) %>%
sf::st_as_sf() %>%
sf::st_cast("POLYGON") %>%
dplyr::mutate(layer = dplyr::row_number()) %>%
dplyr::mutate(
openining_area_m2 = sf::st_area(geometry) %>% as.numeric()
, tree_clump_dist_m = tree_clump_dist_m
)
# return
return(list(dist_rast = dist_rast, openings_vect = openings_vect))
}
# get it
dist_rast_temp = get_clump_dist_rast(clump_polys_temp)what is it?
dist_rast_temp## $dist_rast
## class : SpatRaster
## dimensions : 708, 708, 1 (nrow, ncol, nlyr)
## resolution : 0.2, 0.2 (x, y)
## extent : 341585.1, 341726.7, 4006682, 4006824 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / UTM zone 13N (EPSG:26913)
## source(s) : memory
## name : layer
## min value : 0.00000
## max value : 13.64844
##
## $openings_vect
## Simple feature collection with 52 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 341585.1 ymin: 4006682 xmax: 341726.7 ymax: 4006824
## Projected CRS: NAD83 / UTM zone 13N
## First 10 features:
## layer geometry openining_area_m2 tree_clump_dist_m
## 1 1 POLYGON ((341585.1 4006824,... 0.08 6
## 1.1 2 POLYGON ((341596.5 4006824,... 2.88 6
## 1.2 3 POLYGON ((341660.3 4006824,... 8.32 6
## 1.3 4 POLYGON ((341622.5 4006824,... 7.08 6
## 1.4 5 POLYGON ((341693.1 4006824,... 18.68 6
## 1.5 6 POLYGON ((341676.1 4006818,... 2.48 6
## 1.6 7 POLYGON ((341695.5 4006816,... 0.20 6
## 1.7 8 POLYGON ((341662.5 4006802,... 0.04 6
## 1.8 9 POLYGON ((341645.7 4006804,... 3.40 6
## 1.9 10 POLYGON ((341662.1 4006801,... 0.04 6just view the distance raster
terra::plot(dist_rast_temp$dist_rast)
plot the distance raster and openings vector data we just got with overlaid tree clumps and tree points
plt_fnl_temp =
ggplot() +
# distance
geom_tile(
data = dist_rast_temp$dist_rast %>% terra::aggregate(2, cores = 4) %>% as.data.frame(xy = T) %>% rename(f=3)
, mapping = aes(x=x, y=y, fill = f)
) +
harrypotter::scale_fill_hp(
"mischief"
, na.value = "transparent"
, name = "distance to\nnearest tree (m)"
, limits = c(0,23)
) +
# openings
geom_sf(data = dist_rast_temp$openings_vect, mapping = aes(color = openining_area_m2), fill = NA) +
scale_color_gradient(
low = "gray77", high = "gray11"
, labels = scales::comma_format(accuracy = 1)
, name = latex2exp::TeX("opening\narea ($\\m^2$)")
, limits = c(0,7500)
) +
# clumps
ggnewscale::new_scale_fill() +
geom_sf(
data = clump_polys_temp
, mapping = aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
# tree points
geom_sf(data = ttops_temp, color = "gray88", shape = ".") +
theme_void() +
theme(
plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold")
, legend.title=element_text(size=8), legend.text=element_text(size = 7)
)
# plot
plt_fnl_temp
# save it
ggsave("../data/03_clumps_dist.jpg", width = 7, height = 7)highlight the openings
plt_open_temp =
ggplot() +
# clumps
geom_sf(
data = clump_polys_temp
, mapping = aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
# openings
ggnewscale::new_scale_fill() +
geom_sf(data = dist_rast_temp$openings_vect, mapping = aes(fill = openining_area_m2), color = NA) +
scale_fill_gradient(
low = "gray77", high = "gray11"
, labels = scales::comma_format(accuracy = 1)
, name = latex2exp::TeX("opening\narea ($\\m^2$)")
, limits = c(0,7500)
) +
# tree points
geom_sf(data = ttops_temp, color = "gray88", shape = ".") +
theme_void() +
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
# plot
plt_open_temp
ggsave("../data/04_clumps_openings.jpg", width = 7, height = 7)combine them?
plt_fnl_temp + (plt_open_temp + theme(legend.position = "none")) &
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
ggsave("../data/05_clumps_combined.jpg", width = 10, height = 7)Clump Metrics
create a function to summarize by number of tree clump grouping variable
# create a function to summarize by number of tree clump grouping
get_clump_n_trees_grp_summary = function(trees, clumps){
# get area of harvest unit
#...will use this area in the area calculations such that...
#...TPA = trees in a certain group size across the whole stand area
harvest_area_ha = stand_sf %>%
sf::st_as_sf() %>%
dplyr::mutate(harvest_area_ha = sf::st_area(.) %>% as.numeric() %>% `/`(10000)) %>%
dplyr::pull(harvest_area_ha) %>%
.[1]
# collapse and calculate silv metrics
dta = trees %>%
sf::st_drop_geometry() %>%
dplyr::mutate(stand_area_ha = harvest_area_ha) %>%
dplyr::group_by(stand_area_ha,clump_n_trees_grp) %>%
dplyr::summarise(
# summary metrics
n_trees = dplyr::n_distinct(treeID)
, mean_dbh_cm = mean(dbh_cm, na.rm = T)
, mean_tree_height_m = mean(tree_height_m, na.rm = T)
, loreys_height_m = sum(basal_area_m2*tree_height_m, na.rm = T) / sum(basal_area_m2, na.rm = T)
, basal_area_m2 = sum(basal_area_m2, na.rm = T)
, sum_dbh_cm_sq = sum(dbh_cm^2, na.rm = T)
, .groups = "drop_last"
) %>%
dplyr::ungroup() %>%
# attach clump area
dplyr::left_join(
clumps %>%
sf::st_drop_geometry() %>%
dplyr::group_by(clump_n_trees_grp) %>%
dplyr::summarise(
clump_area_ha = sum(clump_area_ha)
, stand_n_clumps = dplyr::n()
, .groups = "drop_last"
) %>%
dplyr::ungroup()
, by = dplyr::join_by(clump_n_trees_grp)
) %>%
dplyr::mutate(
trees_per_ha = (n_trees/stand_area_ha) # (n_trees/clump_area_ha) ... this was not right
, basal_area_m2_per_ha = (basal_area_m2/stand_area_ha) # (basal_area_m2/clump_area_ha) ... this was not right
, qmd_cm = sqrt(sum_dbh_cm_sq/n_trees)
# stand calcs
, stand_trees_per_ha = sum(n_trees)/stand_area_ha
, stand_basal_area_m2 = sum(basal_area_m2)
, stand_basal_area_m2_per_ha = sum(basal_area_m2)/stand_area_ha
, pct_stand_basal_area = basal_area_m2/stand_basal_area_m2
, pct_stand_n_trees = n_trees/sum(n_trees)
, stand_qmd_cm = sqrt(sum(trees$dbh_cm^2, na.rm = T)/sum(n_trees))
) %>%
dplyr::select(-c(sum_dbh_cm_sq)) %>%
# convert to imperial units
calc_imperial_units_fn()
# return
return(dta)
}call it
# call it
clump_n_trees_grp_summary_temp = get_clump_n_trees_grp_summary(
trees = get_tree_clumps(tree_clump_dist_m = tree_clump_dist_m, ostory_dbh_cm = ostory_dbh_cm)
, clumps = get_clump_summary(
get_tree_clumps(tree_clump_dist_m = tree_clump_dist_m, ostory_dbh_cm = ostory_dbh_cm)
)
)
# what?
clump_n_trees_grp_summary_temp %>% dplyr::glimpse()## Rows: 6
## Columns: 31
## $ stand_area_ha <dbl> 2, 2, 2, 2, 2, 2
## $ clump_n_trees_grp <ord> Individual, 2-4 trees, 5-9 trees, 10-15 tr…
## $ n_trees <int> 12, 30, 31, 11, 40, 645
## $ mean_dbh_cm <dbl> 61.66965, 45.29162, 34.25379, 29.24668, 46…
## $ mean_tree_height_m <dbl> 20.18100, 16.11127, 13.89961, 12.77318, 16…
## $ loreys_height_m <dbl> 24.45865, 24.82122, 20.04323, 15.45925, 23…
## $ basal_area_m2 <dbl> 4.2887219, 7.2452159, 3.8558698, 0.8625523…
## $ clump_area_ha <dbl> 0.03391370, 0.06613346, 0.05884907, 0.0159…
## $ stand_n_clumps <int> 12, 13, 5, 1, 2, 6
## $ trees_per_ha <dbl> 6.0, 15.0, 15.5, 5.5, 20.0, 322.5
## $ basal_area_m2_per_ha <dbl> 2.1443610, 3.6226079, 1.9279349, 0.4312761…
## $ qmd_cm <dbl> 67.45721, 55.45237, 39.79563, 31.59741, 54…
## $ stand_trees_per_ha <dbl> 384.5, 384.5, 384.5, 384.5, 384.5, 384.5
## $ stand_basal_area_m2 <dbl> 101.8253, 101.8253, 101.8253, 101.8253, 10…
## $ stand_basal_area_m2_per_ha <dbl> 50.91266, 50.91266, 50.91266, 50.91266, 50…
## $ pct_stand_basal_area <dbl> 0.042118420, 0.071153376, 0.037867491, 0.0…
## $ pct_stand_n_trees <dbl> 0.01560468, 0.03901170, 0.04031209, 0.0143…
## $ stand_qmd_cm <dbl> 41.06008, 41.06008, 41.06008, 41.06008, 41…
## $ mean_dbh_in <dbl> 24.29784, 17.84490, 13.49599, 11.52319, 18…
## $ qmd_in <dbl> 26.57814, 21.84823, 15.67948, 12.44938, 21…
## $ stand_qmd_in <dbl> 16.17767, 16.17767, 16.17767, 16.17767, 16…
## $ mean_tree_height_ft <dbl> 66.19368, 52.84495, 45.59073, 41.89604, 54…
## $ loreys_height_ft <dbl> 80.22437, 81.41361, 65.74180, 50.70633, 76…
## $ basal_area_ft2_per_ac <dbl> 9.347269, 15.790948, 8.403868, 1.879933, 2…
## $ stand_basal_area_ft2_per_ac <dbl> 221.9283, 221.9283, 221.9283, 221.9283, 22…
## $ trees_per_ac <dbl> 2.4300, 6.0750, 6.2775, 2.2275, 8.1000, 13…
## $ stand_trees_per_ac <dbl> 155.7225, 155.7225, 155.7225, 155.7225, 15…
## $ stand_area_ac <dbl> 4.942, 4.942, 4.942, 4.942, 4.942, 4.942
## $ clump_area_ac <dbl> 0.08380075, 0.16341579, 0.14541606, 0.0393…
## $ basal_area_ft2 <dbl> 46.163803, 77.987504, 41.504582, 9.284513,…
## $ stand_basal_area_ft2 <dbl> 1096.048, 1096.048, 1096.048, 1096.048, 10…summary table
# table it
clump_n_trees_grp_summary_temp %>%
dplyr::select(
clump_n_trees_grp, n_trees
, mean_dbh_in
, qmd_in
, mean_tree_height_ft
, loreys_height_ft
, trees_per_ac
, basal_area_ft2_per_ac, pct_stand_basal_area, pct_stand_n_trees
) %>%
dplyr::mutate(
dplyr::across(
.cols = c(pct_stand_basal_area, pct_stand_n_trees)
, .fns = ~ scales::percent(.x, accuracy = 1)
)
) %>%
kableExtra::kbl(
digits = 1
, escape = F
, caption = paste0("Overstory tree clump summary")
, col.names = c(
"", "trees"
, "mean<br>DBH (in)"
, "QMD (in)"
, "mean<br>Ht. (ft)"
, "Loreys<br>Ht. (ft)"
, "TPA"
, "BA<br>ft<sup>2</sup> ac<sup>-1</sup>"
, "%<br>stand BA"
, "%<br>stand trees"
)
) %>%
kableExtra::kable_styling()| trees |
mean DBH (in) |
QMD (in) |
mean Ht. (ft) |
Loreys Ht. (ft) |
TPA |
BA ft2 ac-1 |
% stand BA |
% stand trees |
|
|---|---|---|---|---|---|---|---|---|---|
| Individual | 12 | 24.3 | 26.6 | 66.2 | 80.2 | 2.4 | 9.3 | 4% | 2% |
| 2-4 trees | 30 | 17.8 | 21.8 | 52.8 | 81.4 | 6.1 | 15.8 | 7% | 4% |
| 5-9 trees | 31 | 13.5 | 15.7 | 45.6 | 65.7 | 6.3 | 8.4 | 4% | 4% |
| 10-15 trees | 11 | 11.5 | 12.4 | 41.9 | 50.7 | 2.2 | 1.9 | 1% | 1% |
| 16-25 trees | 40 | 18.3 | 21.4 | 54.6 | 76.7 | 8.1 | 20.2 | 9% | 5% |
| >25 trees | 645 | 13.6 | 15.3 | 46.2 | 61.6 | 130.6 | 166.3 | 75% | 84% |
ICO Implementation
Churchill et al. (2016) describe the full process for implementing the ICO approach in The ICO Approach to Quantifying and Restoring Forest Spatial Pattern: Implementation Guide in which the authors lay out the prescription development process:
- Identify skips and other special treatment areas
- Consider the need for openings
- Determine the stand average density target
- Determine the appropriate distance to define clumps
- Obtain targets for clump proportions
- Select target clump proportions for your stand
- Generate clump targets for the whole unit
- Combine clump and opening targets with leave tree criteria into marking guidelines
The objective here is to: 1) provide the manager with the current conditions (completed above); 2) take the “targets” as set by the manager (steps 3, 5, 6, 7); 3) create the prescription with the leave tree marking.
Let’s implement this prescription development process with our UAS tree list
3. Determine the stand average density target
Step 3 in Churchill et al. (2016):
An average BA, TPA, or SDI target for the stand should be selected that is appropriate for the species, structure, site conditions, and management objectives. Expected mortality from prescribed fire should be factored in. Stand average targets can come from historical reference stands, plant association based stocking guides, density management tools, or a combination of both (see Franklin et al. (2013) for a full discussion of setting density targets). In dry forests, the number and size of old trees must be accounted in setting the density target. To use the ICO method, the target must be converted to TPA (see Table 1). A lower diameter cutoff also needs to be specified for the TPA target. This should be the lower limit in the contract or cutting guidelines given to the marking crew or contractor. (p.11)
this is what Table 1 looks like with TPA values are derived from the formula:
\[ TPA = \frac{BA}{QMD^{2} \times 0.005454} \]
# function to get tpa from ba and qmd
get_tpa = function(ba_ft2_ac, qmd_in){
tpa = round(ba_ft2_ac/((qmd_in^2)*0.005454))
return(tpa)
}
# table it
tidyr::crossing(
ba = seq(40,200,20)
, qmd = seq(8,20,2)
) %>%
dplyr::mutate(
tpa = get_tpa(ba,qmd)
) %>%
tidyr::pivot_wider(names_from = ba, values_from = tpa) %>%
dplyr::mutate(l = "QMD (in)") %>%
dplyr::relocate(l) %>%
kableExtra::kbl(
col.names = c(".","", seq(40,200,20))
, escape = F
, caption = "Basal Area and QMD to TPA conversion chart"
) %>%
kableExtra::add_header_above(
c("","", "Basal Area (ft2/ac)"=length(seq(40,200,20)))
) %>%
kableExtra::kable_styling() %>%
kableExtra::column_spec(1:2, bold = T) %>%
kableExtra::collapse_rows(columns = 1, valign = "middle")| . | 40 | 60 | 80 | 100 | 120 | 140 | 160 | 180 | 200 | |
|---|---|---|---|---|---|---|---|---|---|---|
| QMD (in) | 8 | 115 | 172 | 229 | 286 | 344 | 401 | 458 | 516 | 573 |
| 10 | 73 | 110 | 147 | 183 | 220 | 257 | 293 | 330 | 367 | |
| 12 | 51 | 76 | 102 | 127 | 153 | 178 | 204 | 229 | 255 | |
| 14 | 37 | 56 | 75 | 94 | 112 | 131 | 150 | 168 | 187 | |
| 16 | 29 | 43 | 57 | 72 | 86 | 100 | 115 | 129 | 143 | |
| 18 | 23 | 34 | 45 | 57 | 68 | 79 | 91 | 102 | 113 | |
| 20 | 18 | 28 | 37 | 46 | 55 | 64 | 73 | 83 | 92 |
Current Stand Conditions
For determining targets, the silviculturist needs to know the current conditions. Provide the current stand conditions based on the UAS tree list for the selected stand that are required to set the targets:
- Current BA
- Current QMD
- Current proportion of trees by clump size
clump_n_trees_grp_summary_temp %>%
dplyr::select(clump_n_trees_grp, pct_stand_n_trees) %>%
dplyr::mutate(
pct_stand_n_trees = scales::percent(pct_stand_n_trees,accuracy = 1)
) %>%
tidyr::pivot_wider(names_from = clump_n_trees_grp, values_from = pct_stand_n_trees) %>%
kableExtra::kable(
caption = paste0(
"Current stand BA (ft2/ac): "
, clump_n_trees_grp_summary_temp$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand QMD (in): "
, clump_n_trees_grp_summary_temp$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand TPA: "
, clump_n_trees_grp_summary_temp$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
)
, escape = F
, digits = 1
) %>%
kableExtra::kable_styling() %>%
kableExtra::footnote(general = "values are the percent of trees in each clump size")| Individual | 2-4 trees | 5-9 trees | 10-15 trees | 16-25 trees | >25 trees |
|---|---|---|---|---|---|
| 2% | 4% | 4% | 1% | 5% | 84% |
| Note: | |||||
| values are the percent of trees in each clump size |
5. Obtain targets for clump proportions
Step 5 in Churchill et al. (2016):
ICO prescriptions are based on a target proportion of trees in different sized clumps within a stand. Proportions are just the percentage of trees, or TPA, that are in different sized clumps. Basal area proportions can be used, but we have found TPA targets to be more straightforward to use. Ideally, a table summarizing clump proportions for a range of reference conditions in your area is available (Table 2). If not, instructions for developing one are provided in section VI. (p.12)
Section VI of Churchill et al. (2016) notes that
reference spatial information may already be available and summarized in a way that it can be directly incorporated into ICO prescriptions. Such data exist and have been published for areas in Arizona (Abella and Denton 2009, Sánchez Meador et al. 2011), the eastern Washington Cascades (Churchill et al. 2013), the northern Rockies (Larson et al. 2012), and the Sierra Nevada (Lydersen et al. 2013). Reference datasets for using ICO in other forest types, such as coastal Douglas-fir or Pacific silver fir, also exist (Larson and Churchill 2008). (p.28)
Table 2 is:
6. Select target clump proportions for your stand
Now set the desired BA, QMD, and proportion of trees in each clump size:
########################################################################################
########################################################################################
# desired BA, QMD, and proportion of trees in each clump size
########################################################################################
########################################################################################
# desired BA
target_ba = 85 # cannot be > current BA
# desired QMD
target_qmd = 18
# desired proportion (%) of trees in each clump size
# !cannot be create larger proportion of ">25 trees" clump as this would require adding trees...
# c("Individual", "2-4 trees", "5-9 trees", "10-15 trees","16-25 trees",">25 trees")
# c(.18, .33, .24, .10, .15)
target_pcts = c(30, 25, 20, 10, 15, 0)
########################################################################################
########################################################################################
# desired BA, QMD, and proportion of trees in each clump size
########################################################################################
########################################################################################Check set up and define data with targets based on Churchill et al. (2016) and their prescription worksheet Excel file to help develop ICO prescriptions
get_target_check_prescription = function(
clump_n_trees_grp_summary_dta
, target_ba = as.numeric(NA)
, target_qmd = as.numeric(NA)
, target_pcts = as.numeric(NA)
){
if(
is.na(target_ba) | is.na(target_qmd) | max(is.na(target_pcts))==1
){
stop("must set all of the function parameters:\n`target_ba`, `target_qmd`, and `target_pcts`")
}
#############################################
# check target BA and TPA
#############################################
if(as.numeric(target_ba)>clump_n_trees_grp_summary_dta$stand_basal_area_ft2_per_ac[1]){
stop(
"target BA in `target_ba` of "
, round(as.numeric(target_ba),1), " is greater than current BA of "
, clump_n_trees_grp_summary_dta$stand_basal_area_ft2_per_ac[1] %>% round(1)
)
}
if(
get_tpa(target_ba, target_qmd)>clump_n_trees_grp_summary_dta$stand_trees_per_ac[1]
){
stop(
"target TPA in of "
, round(as.numeric(get_tpa(target_ba, target_qmd)),1), " is greater than current TPA of "
, clump_n_trees_grp_summary_dta$stand_trees_per_ac[1] %>% round(1)
, "\n adjust `target_ba` and/or `target_qmd` to get valid TPA"
)
}
#############################################
# define data with current and target
# ... this is "smart" in that percentages are adj based on:
# ... 0) are there missing targets?
# ... ... if < 6 numbers provided in `target_pcts` then the largest tree groups get targets of 0
# ... 1) do targets sum to 1?
# ... ... if not trees are distributed proportionally based on targets provided and trees available
# ... 2) is target in largest clump size > current conditions?
# ... ... if yes, target is set to current condition
# ... 3) is target listed in clump size > current largest clump with trees?
# ... ... if yes, target for largest clump size is shifted to current largest clump with trees
#############################################
target_data =
# create data for joining if missing clump groups
dplyr::tibble(
stand_area_ac = rep(clump_n_trees_grp_summary_dta$stand_area_ac[1],6)
, clump_n_trees_grp = factor(
c(1:6)
, labels = c("Individual", "2-4 trees", "5-9 trees", "10-15 trees","16-25 trees",">25 trees")
, ordered = T
)
, mean_clump_n_trees = c(1,3,7,12,20,30)
, min_clump_n_trees = c(1,2,5,10,16,26)
, max_clump_n_trees = c(1,4,9,15,25,99999)
) %>%
dplyr::left_join(
clump_n_trees_grp_summary_dta %>%
dplyr::ungroup() %>%
dplyr::select(clump_n_trees_grp, pct_stand_n_trees, stand_n_clumps)
) %>%
dplyr::mutate(
pct_stand_n_trees = dplyr::coalesce(pct_stand_n_trees,0)
, stand_n_clumps = dplyr::coalesce(stand_n_clumps,0)
) %>%
# add targets
dplyr::bind_cols(
pct_stand_n_trees_target = c(as.numeric(target_pcts), rep(0,6))[1:6] # pad target with 0's
) %>%
# adjust target based on difference from 1
dplyr::mutate(
pct_stand_n_trees_target = pct_stand_n_trees_target*(1/sum(pct_stand_n_trees_target))
# largest clump size with trees
, largest_w_trees = max(ifelse(dplyr::coalesce(pct_stand_n_trees)>0,clump_n_trees_grp,NA),na.rm = T)
, largest_w_trees_target = max(ifelse(dplyr::coalesce(pct_stand_n_trees_target)>0,clump_n_trees_grp,NA),na.rm = T)
) %>%
# move target for largest clump size to the largest current clump size
dplyr::mutate(
pct_stand_n_trees_target = dplyr::case_when(
as.numeric(clump_n_trees_grp)==largest_w_trees &
largest_w_trees_target>largest_w_trees ~ max(
ifelse(as.numeric(clump_n_trees_grp)==largest_w_trees_target,pct_stand_n_trees_target,0)
)
, T ~ pct_stand_n_trees_target
)
) %>%
# adjust target based on current conditions
dplyr::mutate(
pct_stand_n_trees_target = dplyr::case_when(
as.numeric(clump_n_trees_grp)>largest_w_trees &
pct_stand_n_trees_target > 0 ~ 0
, as.numeric(clump_n_trees_grp)==largest_w_trees &
pct_stand_n_trees_target > pct_stand_n_trees ~ pct_stand_n_trees
, T ~ pct_stand_n_trees_target
)
) %>%
# finally, re-scale again based on adjustments
dplyr::mutate(
pct_stand_n_trees_target = dplyr::case_when(
as.numeric(clump_n_trees_grp)==largest_w_trees ~ pct_stand_n_trees_target
, T ~ pct_stand_n_trees_target * (
# pct remaining to scale to
(1-max(ifelse(as.numeric(clump_n_trees_grp)==largest_w_trees,pct_stand_n_trees_target,0))) /
# current pct remaining total allocated
sum(
ifelse(as.numeric(clump_n_trees_grp)!=largest_w_trees,pct_stand_n_trees_target,0))
)
)
) %>%
# add other targets
dplyr::rename(pct_stand_n_trees_current = pct_stand_n_trees) %>%
dplyr::mutate(
stand_trees_per_ac_current = clump_n_trees_grp_summary_dta$stand_trees_per_ac[1]
, stand_trees_per_ac_target = dplyr::coalesce(get_tpa(target_ba, target_qmd),0)
, trees_per_acre_current = stand_trees_per_ac_current*pct_stand_n_trees_current
, trees_per_acre_target = stand_trees_per_ac_target*pct_stand_n_trees_target
, clumps_per_acre_current = trees_per_acre_current/mean_clump_n_trees
, clumps_per_acre_target = trees_per_acre_target/mean_clump_n_trees
, stand_n_clumps_current = stand_n_clumps
, stand_n_clumps_target = (clumps_per_acre_target*stand_area_ac) %>% round(0)
) %>%
dplyr::select(-c(tidyselect::starts_with("largest_w_trees"), stand_n_clumps))
# ????
# target_data %>% glimpse()
# issue warning about targets
if(min(target_data$pct_stand_n_trees_target == c(as.numeric(target_pcts), rep(0,6))[1:6])==0){
warning(
"proportion of trees in each clump size target `target_pcts` adjusted!!!"
, "\nfrom : ", paste(round(target_pcts,2),collapse = ",")
, "\nto : ", paste(round(target_data$pct_stand_n_trees_target,2),collapse = ",")
)
}
# return
return(target_data)
}
# call it
target_data_temp = get_target_check_prescription(
clump_n_trees_grp_summary_temp
, target_ba = target_ba
, target_qmd = target_qmd
, target_pcts = target_pcts
)## Warning in get_target_check_prescription(clump_n_trees_grp_summary_temp, : proportion of trees in each clump size target `target_pcts` adjusted!!!
## from : 30,25,20,10,15,0
## to : 0.3,0.25,0.2,0.1,0.15,0# what?
target_data_temp %>% dplyr::glimpse()## Rows: 6
## Columns: 15
## $ stand_area_ac <dbl> 4.942, 4.942, 4.942, 4.942, 4.942, 4.942
## $ clump_n_trees_grp <ord> Individual, 2-4 trees, 5-9 trees, 10-15 tre…
## $ mean_clump_n_trees <dbl> 1, 3, 7, 12, 20, 30
## $ min_clump_n_trees <dbl> 1, 2, 5, 10, 16, 26
## $ max_clump_n_trees <dbl> 1, 4, 9, 15, 25, 99999
## $ pct_stand_n_trees_current <dbl> 0.01560468, 0.03901170, 0.04031209, 0.01430…
## $ pct_stand_n_trees_target <dbl> 0.30, 0.25, 0.20, 0.10, 0.15, 0.00
## $ stand_trees_per_ac_current <dbl> 155.7225, 155.7225, 155.7225, 155.7225, 155…
## $ stand_trees_per_ac_target <dbl> 48, 48, 48, 48, 48, 48
## $ trees_per_acre_current <dbl> 2.4300, 6.0750, 6.2775, 2.2275, 8.1000, 130…
## $ trees_per_acre_target <dbl> 14.4, 12.0, 9.6, 4.8, 7.2, 0.0
## $ clumps_per_acre_current <dbl> 2.4300000, 2.0250000, 0.8967857, 0.1856250,…
## $ clumps_per_acre_target <dbl> 14.400000, 4.000000, 1.371429, 0.400000, 0.…
## $ stand_n_clumps_current <dbl> 12, 13, 5, 1, 2, 6
## $ stand_n_clumps_target <dbl> 71, 20, 7, 2, 2, 0current vs target
target_data_temp %>%
dplyr::select(
clump_n_trees_grp
, tidyselect::starts_with("pct_stand_n_trees")
, tidyselect::starts_with("trees_per_acre_")
) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("pct_stand_n_trees")
, ~ scales::percent(.x,accuracy = 1)
)
) %>%
kableExtra::kable(
caption = paste0(
"Current stand BA (ft2/ac): "
, clump_n_trees_grp_summary_temp$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand QMD (in): "
, clump_n_trees_grp_summary_temp$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand TPA: "
, clump_n_trees_grp_summary_temp$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
)
, escape = F
, digits = 1
, col.names = c(
"", rep(c("current","target"),2)
)
) %>%
kableExtra::kable_styling() %>%
kableExtra::add_header_above(
c(" "=1,"% Trees"=2, "TPA"=2)
)| current | target | current | target | |
|---|---|---|---|---|
| Individual | 2% | 30% | 2.4 | 14.4 |
| 2-4 trees | 4% | 25% | 6.1 | 12.0 |
| 5-9 trees | 4% | 20% | 6.3 | 9.6 |
| 10-15 trees | 1% | 10% | 2.2 | 4.8 |
| 16-25 trees | 5% | 15% | 8.1 | 7.2 |
| >25 trees | 84% | 0% | 130.6 | 0.0 |
8. Combine clump and opening targets with leave tree criteria into marking guidelines
Use our UAS tree list to generate the prescription:
- start with the largest clump size currently with trees
- cut trees to the next largest clump size until desired # clumps is reached
- repeat with each successive clump size through to individual tree selection
- if possible, cut in same clump until desired proportions are reached to minimize machine time
Function to cut clumps
put the entire process outlined immediately above into a function to use with the tree list to cut each clump and then we’ll select the desired proportion of clumps
define intermediate functions
define intermediate functions
##############################################
# working with sf LINESTRINGS
##############################################
# first two functions borrowed from https://github.com/metafor-ulaval/ALSroads/blob/main/R/line_tools.R
########
# Get heading of both ends of a line
########
st_ends_heading <- function(line){
M <- sf::st_coordinates(line)
i <- c(2, nrow(M) - 1)
j <- c(1, -1)
headings <- mapply(i, j, FUN = function(i, j) {
Ax = M[i-j,1]
Ay = M[i-j,2]
Bx = M[i,1]
By = M[i,2]
atan2(Ay-By, Ax-Bx)*180/pi
})
names(headings) <- c("head", "tail")
return(headings)
}
########
# extend the line on both ends
########
st_extend_line <- function(line, distance, end = "BOTH"){
if (!(end %in% c("BOTH", "HEAD", "TAIL")) | length(end) != 1) stop("'end' must be 'BOTH', 'HEAD' or 'TAIL'")
M <- sf::st_coordinates(line)[,-3]
keep <- !(end == c("TAIL", "HEAD"))
ends <- c(1, nrow(M))[keep]
headings <- st_ends_heading(line)[keep] / 180 * pi
distances <- if (length(distance) == 1) rep(distance, 2) else distance[1:2]
M[ends, 1:2] <- M[ends, 1:2] + distances[keep] * c(cos(headings), sin(headings))
newline <- sf::st_linestring(M)
# If input is sfc_LINESTRING and not sfg_LINESTRING
if (is.list(line)) newline <- sf::st_sfc(newline, crs = sf::st_crs(line))
return(newline)
}
########
# pass an sf dataframe of points and return a line between the farthest points
########
st_points_to_line <- function(pts, line_ext=0) {
if(max(class(ttops_temp) %in% c("sf"))!=1){
stop("must provide an object of class `sf`")
}
# find farthest distance between points
dist_temp = sf::st_distance(pts)
# get the points
f_pts_temp =
pts %>%
dplyr::ungroup() %>%
dplyr::slice(
# get the farthest points from distance matrix
which(dist_temp == max(dist_temp), arr.ind = TRUE)[1,]
)
# draw a line between the farthest two points
f_line_temp = f_pts_temp %>%
# convert to linestring
dplyr::ungroup() %>%
dplyr::summarise(n=dplyr::n()) %>%
sf::st_cast("LINESTRING") %>%
dplyr::ungroup() %>%
dplyr::select(-c(n))
# and apply the line extension
farthest_line = st_extend_line(f_line_temp, distance = line_ext)
# return
return(farthest_line)
}
# st_points_to_line(ttops_temp %>% dplyr::slice_head(prop = .1), line_ext = 6) %>%
# ggplot() + geom_sf() + geom_sf(data = ttops_temp %>% dplyr::slice_head(prop = .1)) + theme_void()
########
# Function to calculate Euclidean distance between 2 points
########
st_euclidean_distance <- function(p1,p2) {
return(sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2))
}
########
# return a line perpendicular to current line
########
### https://stackoverflow.com/questions/56771058/perpendicular-lines-at-regular-intervals-along-lines-with-multiple-nodes
# Function to calculate 2 points on a line perpendicular to another defined by 2 points p1,p2
# For point at interval, which can be a proportion of the segment length, or a constant
st_perp_line <- function(interval=0.5, my_line, proportion=TRUE) {
# get end points of line
p1 = my_line %>% sf::st_cast("POINT") %>% sf::st_coordinates() %>% .[1,]
p2 = my_line %>% sf::st_cast("POINT") %>% sf::st_coordinates() %>% .[2,]
# get length of line to return equal length line
line_len = sf::st_length(my_line) %>% as.numeric() %>% `/`(2)
# get crs of line
my_crs = sf::st_crs(my_line)
# Calculate x and y distances
x_len <- p2[1] - p1[1]
y_len <- p2[2] - p1[2]
# If proportion calculate reference point from tot_length
if (proportion) {
point <- c(p1[1]+x_len*interval,p1[2]+y_len*interval)
}
# Else use the constant value
else {
tot_len <- st_euclidean_distance(p1,p2)
point <- c(p1[1]+x_len/tot_len*interval,p1[2]+y_len/tot_len*interval)
}
# Calculate the x and y distances from reference point to point on line line_len distance away
ref_len <- st_euclidean_distance(point,p2)
xn_len <- (line_len / ref_len) * (p2[1] - point[1])
yn_len <- (line_len / ref_len) * (p2[2] - point[2])
# fix for identical
if(identical(point,p2) & x_len>y_len){ # this works for horizontal line
xn_len <- line_len/2
yn_len <- 0
}else if(identical(point,p2) & x_len<y_len){ # this works for vertical line
xn_len <- 0
yn_len <- line_len/2
}
# Invert the x and y lengths and add/subtract from the refrence point
# ref_points <- rbind(point,c(point[1] + yn_len,point[2] - xn_len),c(point[1] - yn_len,point[2] + xn_len))
ref_points <- rbind(c(point[1] + yn_len,point[2] - xn_len),c(point[1] - yn_len,point[2] + xn_len))
# use the reference points to return a line
return(
ref_points %>%
dplyr::as_tibble() %>%
dplyr::rename_with(tolower) %>%
sf::st_as_sf(coords = c("x", "y"), crs = my_crs, remove = F) %>%
dplyr::summarise(n=dplyr::n()) %>%
sf::st_cast("LINESTRING") %>%
dplyr::ungroup() %>%
dplyr::select(-c(n))
)
}function to pass a clump
generate data using clumping functions above, pass that data, return tree list with cut/keep flag based on the target clump group size
# function to cut a single clump
cut_clump_fn <- function(
c # clump id from need_cut_trees data
, need_cut_trees # data with clumps already defined returned by st_clump_points()
, tgt = tgt # select one level of input passed to st_clump_points() ...
# ... clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees")
){
# use only the current clump
curr_dta = need_cut_trees %>%
dplyr::filter(clump_id == c)
if(nrow(curr_dta)<1){
stop("cannot find data with the clump_id == `c` in `need_cut_trees` data")
}
# get distance used to create clumps in st_clump_points()
dist_temp = curr_dta$tree_clump_dist_m[1]
# do we even need to cut?
reclump = st_clump_points(x = curr_dta, clump_dist_m = dist_temp)
if(
(
reclump %>%
dplyr::pull(clump_id) %>%
unique() %>%
length()
) != 1
){
stop("this is not a clump...send the `need_cut_trees` data through the st_clump_points function again")
}
# do we even need to cut?
if(
unique(reclump$clump_n_trees_grp) == tgt
){
return(
reclump %>%
sf::st_drop_geometry() %>%
dplyr::select(treeID, is_keep_tree)
)
}
# create clump polygon
clumps = get_clump_summary(curr_dta)
# get the farthest line between points
f_line_temp = st_points_to_line(curr_dta, line_ext = dist_temp)
# get perpendicular lines in dataset which we can iterate over to make cuts
perp_line_sf_temp =
# for every 1 m along line length, get a new perp line
seq(
from = 0
, to = sf::st_length(f_line_temp) %>% as.numeric() %>% floor()
, by = 1
) %>%
purrr::map(
st_perp_line
, my_line = f_line_temp
, proportion = F
) %>%
dplyr::bind_rows() %>%
dplyr::mutate(line_n = dplyr::row_number())
# find intersection of lines with the polygon and add length of intersection to perp line data
perp_line_sf_temp = perp_line_sf_temp %>%
dplyr::inner_join(
# intersect and calc len
perp_line_sf_temp %>%
sf::st_intersection(
clumps %>%
dplyr::ungroup() %>%
dplyr::select(clump_id) %>%
dplyr::filter(clump_id == c)
) %>%
dplyr::mutate(len_m = sf::st_length(geometry) %>% as.numeric()) %>%
sf::st_drop_geometry()
, by = "line_n"
)
# make cuts at the points where there is the least overlap with the clump polygon
# list of potential line cut + tree combinations
# aggregate the mean length of the intersecting cut lines to the tree level
# prioritize trees for removal that have smallest length
cut_tree_lines_temp =
curr_dta %>%
sf::st_buffer(dist_temp/2) %>%
sf::st_join(perp_line_sf_temp) %>%
sf::st_drop_geometry() %>%
dplyr::group_by(treeID) %>%
dplyr::summarise(len_m = mean(len_m, na.rm = T)) %>%
dplyr::ungroup() %>%
dplyr::arrange(len_m, treeID) %>%
dplyr::mutate(n = dplyr::row_number())
###############################################################
# while
###############################################################
# cut until the desired clump size is achieved based on these cut lines
while_temp = 1
i = 1
while(while_temp==1) {
# cut trees
cut_trees_temp = cut_tree_lines_temp %>%
dplyr::slice(1:i) %>%
dplyr::distinct(treeID)
# get the remaining trees not cut
trees_remain_temp = curr_dta %>%
dplyr::anti_join(cut_trees_temp, by = "treeID")
# ensure that there are trees
if(nrow(trees_remain_temp)==0){
if(best_desired_grps_n_temp==0){
# there are no more possible cuts :/
# ... so we're going to add trees until the desired clump size is reached
# start with biggest tree and keep adding trees until the desired clump size is reached
start_tree = curr_dta %>%
dplyr::arrange(desc(dbh_cm)) %>%
dplyr::filter(dplyr::row_number() == 1)
increment_m = 0.5
k = 1
keep_trees = dplyr::tibble(treeID = character(0))
while(nrow(keep_trees)==0){
# create a polygon to intersect with tree points and keep trees until number of trees met
poly_keep = start_tree %>%
sf::st_buffer( (dist_temp/2) + increment_m*k ) %>%
dplyr::mutate(dummy = 1) %>%
dplyr::select(dummy)
# intersect and clump
i_trees = curr_dta %>%
sf::st_intersection(poly_keep) %>%
st_clump_points(clump_dist_m = dist_temp)
i_trees %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp, clump_id)
# keep only the desired clump
keep_trees = i_trees %>%
sf::st_drop_geometry() %>%
dplyr::filter(
clump_id ==
i_trees %>%
dplyr::filter(clump_n_trees_grp == tgt) %>%
dplyr::pull(clump_id) %>%
.[1] %>%
dplyr::coalesce("nope")
) %>%
dplyr::select(treeID)
# increment
k = k+1
}
# the remaining trees are cuts
best_cuts = curr_dta %>%
sf::st_drop_geometry() %>%
dplyr::anti_join(keep_trees, by = "treeID") %>%
dplyr::select(treeID)
}else if(
best_desired_grps_n_temp>0
){
###############################
# add trees back in until gets worse
###############################
j = 1
while_add = 1
while(while_add==1){
# cut trees
cut_trees_temp = best_cuts %>%
# add trees back in (i.e. remove them from the cut trees)
dplyr::anti_join(
cut_tree_lines_temp %>%
dplyr::slice(1:j) %>%
dplyr::distinct(treeID)
, by = "treeID"
)
# get the remaining trees not cut
trees_remain_temp = curr_dta %>%
dplyr::anti_join(cut_trees_temp, by = "treeID")
# count the groups remaining after cuts
desired_grps_n_temp = trees_remain_temp %>%
st_clump_points(clump_dist_m = dist_temp) %>%
sf::st_drop_geometry() %>%
# do we have group sizes we want?
dplyr::filter(clump_n_trees_grp == tgt) %>%
nrow()
if(desired_grps_n_temp>=best_desired_grps_n_temp){ # is this better than the best
best_cuts = cut_trees_temp
best_desired_grps_n_temp = desired_grps_n_temp
}else{
# stop it
while_add = 0
}
# increment
j = j + 1
} # while(while_add==1)
} # best_desired_grps_n_temp>0
# done so stop the whole stop it
while_temp = 0
}else{ # if(nrow(trees_remain_temp)==0)
# count the groups remaining after cuts
desired_grps_n_temp = trees_remain_temp %>%
st_clump_points(clump_dist_m = dist_temp) %>%
sf::st_drop_geometry() %>%
# do we have group sizes we want?
dplyr::filter(clump_n_trees_grp == tgt) %>%
nrow()
### store best cut list
if(i==1){
best_cuts = cut_trees_temp
best_desired_grps_n_temp = desired_grps_n_temp
}else if(desired_grps_n_temp>best_desired_grps_n_temp){ # is this better than the best
best_cuts = cut_trees_temp
best_desired_grps_n_temp = desired_grps_n_temp
}else if(
desired_grps_n_temp==best_desired_grps_n_temp
& i!=nrow(cut_tree_lines_temp)
){
best_cuts = best_cuts
best_desired_grps_n_temp = best_desired_grps_n_temp
}else if(
best_desired_grps_n_temp>0
& desired_grps_n_temp<best_desired_grps_n_temp
){ # is this worse than the best which was successful?
###############################
# add trees back in until gets worse
###############################
j = 1
while_add = 1
while(while_add==1){
# cut trees
cut_trees_temp = best_cuts %>%
# add trees back in (i.e. remove them from the cut trees)
dplyr::anti_join(
cut_tree_lines_temp %>%
dplyr::slice(1:j) %>%
dplyr::distinct(treeID)
, by = "treeID"
)
# get the remaining trees not cut
trees_remain_temp = curr_dta %>%
dplyr::anti_join(cut_trees_temp, by = "treeID")
# count the groups remaining after cuts
desired_grps_n_temp = trees_remain_temp %>%
st_clump_points(clump_dist_m = dist_temp) %>%
sf::st_drop_geometry() %>%
# do we have group sizes we want?
dplyr::filter(clump_n_trees_grp == tgt) %>%
nrow()
if(desired_grps_n_temp>=best_desired_grps_n_temp){ # is this better than the best
best_cuts = cut_trees_temp
best_desired_grps_n_temp = desired_grps_n_temp
}else{
# stop it
while_add = 0
}
# increment
j = j + 1
} # while(while_add==1)
# done so stop the whole stop it
while_temp = 0
}else if( i==nrow(cut_tree_lines_temp) ){ # is this the end?
# stop it
while_temp = 0
}
### increment
i = i+1
} # else if(nrow(trees_remain_temp)==0)
} # while(while_temp==1)
# return it
# return treelist with cut/keep
# join to original data and pull
d_temp = curr_dta %>%
sf::st_drop_geometry() %>%
dplyr::mutate(is_keep_tree = as.numeric(NA)) %>%
dplyr::select(-c(is_keep_tree)) %>%
dplyr::left_join(
best_cuts %>% dplyr::mutate(is_keep_tree = 0)
, by = dplyr::join_by("treeID")
) %>%
dplyr::mutate(is_keep_tree = dplyr::coalesce(is_keep_tree, 1)) %>%
dplyr::select(treeID, is_keep_tree)
# returns treeID and keep tree flag data frame
return(d_temp)
}
# # example
# # pass it a clump id in the data
# ttops_temp %>%
# dplyr::filter(clump_n_trees_grp == "2-4 trees") %>%
# dplyr::arrange(desc(clump_n_trees)) %>%
# dplyr::slice(1) %>%
# dplyr::pull(clump_id) %>%
# cut_clump_fn(
# need_cut_trees = ttops_temp
# , tgt = "Individual"
# )there are some cases where the cut_clumps_fn may not
return the desired tree clumps based on a single pass using the initial
cut lines…iterate over the function until all clumps are in the desired
clump size or smaller.
get_cut_clump <- function(
c # clump id from need_cut_trees data
, need_cut_trees # data with clumps already defined returned by st_clump_points()
, tgt = tgt # select one level of input passed to st_clump_points() ...
# ... clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees")
) {
# what if we need to keep cutting because could not find a solution based on initial cut lines?
cuts_first = cut_clump_fn(
c = c
, need_cut_trees = need_cut_trees
, tgt = tgt
)
# did we get the target?
new_clumps =
need_cut_trees %>%
dplyr::inner_join(
cuts_first %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
, by = "treeID"
) %>%
st_clump_points()
# new_clumps %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp)
while(
(
new_clumps %>% dplyr::filter(clump_n_trees_grp > tgt) %>% nrow()
) > 0
){
# redo the cut clump
# redo the cut clump
cuts_again = new_clumps %>%
dplyr::filter(clump_n_trees_grp > tgt) %>%
dplyr::pull(clump_id) %>%
unique() %>%
purrr::map(cut_clump_fn, need_cut_trees = new_clumps, tgt = tgt) %>%
dplyr::bind_rows() %>%
dplyr::rename(updt = is_keep_tree)
# update the original cut data
cuts_first = cuts_first %>%
dplyr::left_join(
cuts_again
, by = "treeID"
) %>%
dplyr::mutate(is_keep_tree = dplyr::coalesce(updt, is_keep_tree)) %>%
dplyr::select(-c(updt))
# reset the new clumps
new_clumps =
need_cut_trees %>%
dplyr::inner_join(
cuts_first %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
, by = "treeID"
) %>%
st_clump_points()
# new_clumps %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp)
} # end while
return(cuts_first)
}function to pass a tree list
function to pass a whole tree list of sf point data
##############################################
# get cut keep tree flag
##############################################
get_keep_tree_flag <- function(
x # x = sf point data
, tgt # tgt = clump_n_trees_grp level defined in st_clump_points: ...
## ... clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees")
, clump_dist_m = 6 # size (radius) of the epsilon neighborhood = maximum distance between points to add to clump
) {
# MAKE CUTS INDEPENDENT OF CURRENT CLUMP GROUPING...
# what if we pass a tree list with trees already cut? and thus, potentially different clump group sizes than is
# currently defined in clump_n_trees_grp?
# 1) apply the clump grouping
# 2) only go through the cut alg for clumps that need cutting to the target
# 3) append the trees in target clump or lower to the list at the end
#############################
# 1) apply the clump grouping
#############################
new_clump_trees = x %>%
st_clump_points(clump_dist_m = clump_dist_m)
# need cutting still
need_cut_trees = new_clump_trees %>%
dplyr::filter(clump_n_trees_grp > tgt)
#############################
# 2) only go through the cut alg for clumps that need cutting to the target
#############################
# for each clump that still needs cutting, iterate over and return keep/cut flag
compl_need_cut_trees = need_cut_trees %>%
dplyr::pull(clump_id) %>%
unique() %>%
purrr::map(get_cut_clump, need_cut_trees = need_cut_trees, tgt = tgt) %>% # purrr:map fn
dplyr::bind_rows()
#############################
# 4) append the trees in target clump or lower to the list at the end
#############################
compl_need_cut_trees = compl_need_cut_trees %>%
dplyr::bind_rows(
new_clump_trees %>%
sf::st_drop_geometry() %>%
dplyr::filter(clump_n_trees_grp <= tgt) %>%
dplyr::mutate(is_keep_tree = 1) %>%
dplyr::select(treeID, is_keep_tree)
)
# return
return(
# original data so that order is preserved
x %>%
# add on trees that got cut with trees already good
dplyr::inner_join(
compl_need_cut_trees
, by = "treeID"
) %>%
dplyr::pull(is_keep_tree)
)
} # get_keep_tree_flagProcess to cut groups incrementally
cut groups incrementally to achieve target tree clump group size proportions
# must first have target data
# target_data_temp = get_target_check_prescription(
# clump_n_trees_grp_summary_temp
# , target_ba = target_ba
# , target_qmd = target_qmd
# , target_pcts = target_pcts
# )
# what is the smallest group with a target?
sm_grp_temp = target_data_temp %>%
dplyr::filter(pct_stand_n_trees_target>0) %>%
dplyr::pull(clump_n_trees_grp) %>%
min()
# start with the largest group in target data currently with trees
grp_temp = target_data_temp %>%
dplyr::arrange(desc(clump_n_trees_grp)) %>%
dplyr::filter(pct_stand_n_trees_current>0) %>%
dplyr::slice(1) %>%
dplyr::pull(clump_n_trees_grp)
message(
paste("started cutting", grp_temp, "at", Sys.time())
)
# identify clumps to leave untouched
# clump_polys_temp = get_clump_summary(ttops_temp)
keep_clumps_temp =
clump_polys_temp %>%
sf::st_drop_geometry() %>%
#keep only trees in current group
dplyr::inner_join(
target_data_temp %>%
dplyr::filter(clump_n_trees_grp == grp_temp) %>%
dplyr::select(
clump_n_trees_grp, mean_clump_n_trees, min_clump_n_trees, max_clump_n_trees
, stand_n_clumps_target
)
, by = "clump_n_trees_grp"
) %>%
# keep only clumps that meet criteria
dplyr::filter(
n_trees >= min_clump_n_trees
& n_trees <= max_clump_n_trees
) %>%
# keep clumps closest the mean number in target
dplyr::mutate(
pct_to_target = abs( (n_trees-mean_clump_n_trees)/mean_clump_n_trees )
) %>%
dplyr::arrange(pct_to_target, desc(qmd_cm), n_trees, clump_id) %>%
dplyr::filter(dplyr::row_number()<=stand_n_clumps_target) %>%
dplyr::select(clump_id)
# !!!!!!!!!!FIX: WHAT IF WE HAVE CLUMPS OF THIS SIZE BUT THE N_TREES>MAX_TREES AND NEED TO GET CLUMPS OF THIS SIZE?
# ... UPDATE CUTTING ALG TO CUT CLUMP DOWN TO MIN-MAX TREE RANGE FOR CLUMPS WITH INF UPPER LIMIT
# start building tree list with keep/cut flag
# ttops_temp = get_tree_clumps(
# my_suid = my_suid
# , tree_clump_dist_m = tree_clump_dist_m
# , ostory_dbh_cm = ostory_dbh_cm
# )
# build tree list
keep_trees =
ttops_temp %>%
dplyr::ungroup() %>%
sf::st_drop_geometry() %>%
dplyr::inner_join(keep_clumps_temp, by = "clump_id") %>%
dplyr::select(treeID) %>%
dplyr::mutate(is_keep_tree = 1)
# did we keep all of the clumps?
if(
nrow(keep_clumps_temp) ==
( clump_polys_temp %>% dplyr::filter(clump_n_trees_grp == grp_temp) %>% nrow() )
){
remaining_trees = dplyr::tibble(
treeID = character(0)
, is_keep_tree = numeric(0)
)
}else{
# determine keep/cut for the remaining trees in that group type
remaining_trees =
ttops_temp %>%
dplyr::ungroup() %>%
# keep only trees in current grp
dplyr::filter(clump_n_trees_grp == grp_temp) %>%
# remove keep trees
dplyr::anti_join(keep_trees, by = "treeID")
# apply the cutting algorithm
# get the flag
remaining_trees$is_keep_tree = get_keep_tree_flag(
x = remaining_trees
, tgt = target_data_temp %>%
dplyr::ungroup() %>%
dplyr::arrange(clump_n_trees_grp) %>%
dplyr::mutate(l = dplyr::lag(clump_n_trees_grp)) %>%
dplyr::filter(clump_n_trees_grp == grp_temp) %>%
dplyr::pull(l)
, clump_dist_m = tree_clump_dist_m
)
# select relevant columns
remaining_trees = remaining_trees %>%
sf::st_drop_geometry() %>%
dplyr::ungroup() %>%
dplyr::select(treeID, is_keep_tree)
}
message(
paste("done cutting", grp_temp, "at", Sys.time())
)
###############################################
# now process to go on to the next groups
###############################################
# get the next group
grp_temp = target_data_temp %>%
dplyr::filter(clump_n_trees_grp<grp_temp) %>%
dplyr::arrange(desc(clump_n_trees_grp)) %>%
dplyr::slice(1) %>%
dplyr::pull(clump_n_trees_grp)
while(grp_temp>=sm_grp_temp & grp_temp != "Individual"){
message(
paste("started cutting", grp_temp, "at", Sys.time())
)
# identify clumps to leave untouched
# clump_polys_temp = get_clump_summary(ttops_temp)
keep_clumps_temp =
clump_polys_temp %>%
sf::st_drop_geometry() %>%
#keep only trees in current group
dplyr::inner_join(
target_data_temp %>%
dplyr::filter(clump_n_trees_grp == grp_temp) %>%
dplyr::select(
clump_n_trees_grp, mean_clump_n_trees, min_clump_n_trees, max_clump_n_trees
, stand_n_clumps_target
)
, by = "clump_n_trees_grp"
) %>%
# keep only clumps that meet criteria
dplyr::filter(
n_trees >= min_clump_n_trees
& n_trees <= max_clump_n_trees
) %>%
# keep clumps closest the mean number in target group
dplyr::mutate(
pct_to_target = abs( (n_trees-mean_clump_n_trees)/mean_clump_n_trees )
) %>%
dplyr::arrange(pct_to_target, desc(qmd_cm), n_trees, clump_id) %>%
dplyr::filter(dplyr::row_number()<=stand_n_clumps_target) %>%
dplyr::select(clump_id)
# add to tree list with keep/cut flag
keep_trees = keep_trees %>%
dplyr::bind_rows(
ttops_temp %>%
dplyr::ungroup() %>%
sf::st_drop_geometry() %>%
dplyr::inner_join(keep_clumps_temp, by = "clump_id") %>%
dplyr::select(treeID) %>%
dplyr::mutate(is_keep_tree = 1)
)
# check if the desired clump number was met and add the remaining trees from previous group if needed
more_clumps_target = 0
if(
nrow(keep_clumps_temp) <
(
target_data_temp %>%
dplyr::filter(clump_n_trees_grp == grp_temp) %>%
dplyr::pull(stand_n_clumps_target)
)
){
# how many more clumps are needed?
more_clumps_target =
(
target_data_temp %>%
dplyr::filter(clump_n_trees_grp == grp_temp) %>%
dplyr::pull(stand_n_clumps_target)
) - nrow(keep_clumps_temp)
# determine group size of remaining trees in the group prior that were not in a group selected and were not cut
potential_trees =
ttops_temp %>%
dplyr::ungroup() %>%
dplyr::inner_join(
remaining_trees %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
, by = "treeID"
) %>%
st_clump_points(clump_dist_m = tree_clump_dist_m) %>%
# ggplot() + geom_sf(aes(fill = clump_n_trees_grp)) + theme_void()
# get the original clump id to prioritize new clump groups in the same area
dplyr::inner_join(
ttops_temp %>%
sf::st_drop_geometry() %>%
dplyr::ungroup() %>%
dplyr::select(treeID, clump_id) %>%
dplyr::rename(orig_clump_id = clump_id)
, by = "treeID"
) %>%
dplyr::group_by(orig_clump_id) %>%
dplyr::mutate(
pct_desired_grp =
sum(ifelse(clump_n_trees_grp == grp_temp, 1, 0)) / dplyr::n()
) %>%
dplyr::ungroup() %>%
# keep only the current group
dplyr::filter(clump_n_trees_grp == grp_temp)
# pick trees from potential trees based on clump summary
keep_trees = keep_trees %>%
dplyr::bind_rows(
potential_trees %>%
sf::st_drop_geometry() %>%
# filter trees based on clumps needed
dplyr::inner_join(
get_clump_summary(potential_trees) %>%
sf::st_drop_geometry() %>%
# join with original clump id metrics to prioritize
# keeping clumps in the same area and minimize cutting time
dplyr::left_join(
potential_trees %>%
sf::st_drop_geometry() %>%
dplyr::group_by(clump_id, orig_clump_id) %>%
dplyr::summarise(pct_desired_grp = max(pct_desired_grp))
, by = "clump_id"
) %>%
# keep the number of clumps needed
dplyr::arrange(desc(pct_desired_grp), orig_clump_id, desc(qmd_cm), n_trees, clump_id) %>%
dplyr::filter(dplyr::row_number()<=more_clumps_target) %>%
dplyr::select(clump_id)
, by = "clump_id"
) %>%
dplyr::select(treeID) %>%
dplyr::mutate(is_keep_tree = 1)
)
} # end if don't have enough clumps
####################################################
# determine keep/cut for the remaining trees in that group type and higher groups
####################################################
remaining_trees =
ttops_temp %>%
dplyr::ungroup() %>%
# REMOVE TREES FROM PREVIOUS TREES REMAINING that got cut to make this group size
dplyr::anti_join(
remaining_trees %>% dplyr::filter(is_keep_tree == 0)
, by = "treeID"
) %>%
# keep only trees in current grp or prior grp
dplyr::filter(clump_n_trees_grp >= grp_temp) %>%
# remove keep trees
dplyr::anti_join(keep_trees, by = "treeID")
if(nrow(remaining_trees) > 0){
# apply the cutting algorithm
# get the flag
remaining_trees$is_keep_tree = get_keep_tree_flag(
x = remaining_trees
, tgt = target_data_temp %>%
dplyr::ungroup() %>%
dplyr::arrange(clump_n_trees_grp) %>%
dplyr::mutate(l = dplyr::lag(clump_n_trees_grp)) %>%
dplyr::filter(clump_n_trees_grp == grp_temp) %>%
dplyr::pull(l)
, clump_dist_m = tree_clump_dist_m
)
# select relevant columns
remaining_trees = remaining_trees %>%
sf::st_drop_geometry() %>%
dplyr::ungroup() %>%
dplyr::select(treeID, is_keep_tree)
}
# increment
message(
paste("done cutting", grp_temp, "at", Sys.time())
)
# get the next group
grp_temp = target_data_temp %>%
dplyr::filter(clump_n_trees_grp<grp_temp) %>%
dplyr::arrange(desc(clump_n_trees_grp)) %>%
dplyr::slice(1) %>%
dplyr::pull(clump_n_trees_grp) %>%
dplyr::coalesce("Individual")
}
# now individuals
if(sm_grp_temp == "Individual"){
message(
paste("started cutting", sm_grp_temp, "at", Sys.time())
)
potential_trees =
# original data
ttops_temp %>%
dplyr::filter(clump_n_trees_grp == "Individual") %>%
# remaining trees
dplyr::bind_rows(
ttops_temp %>%
dplyr::ungroup() %>%
dplyr::inner_join(
remaining_trees %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
, by = "treeID"
)
) %>%
# make sure that these are all individuals
dplyr::mutate(is_orig = ifelse(clump_n_trees_grp=="Individual", 1, 0)) %>%
st_clump_points(clump_dist_m = tree_clump_dist_m) %>%
dplyr::filter(clump_n_trees_grp == "Individual")
# pick trees from potential trees based on target
keep_trees = keep_trees %>%
dplyr::select(treeID) %>%
dplyr::bind_rows(
potential_trees %>%
sf::st_drop_geometry() %>%
# keep the number of clumps needed
dplyr::arrange(desc(is_orig), desc(dbh_cm), desc(tree_height_m)) %>%
dplyr::filter(
dplyr::row_number() <=
target_data_temp %>%
dplyr::filter(clump_n_trees_grp == "Individual") %>%
dplyr::pull(stand_n_clumps_target)
) %>%
dplyr::select(treeID)
) %>%
dplyr::mutate(is_keep_tree = 1)
message(
paste("done cutting", sm_grp_temp, "at", Sys.time())
)
}
# get the final prescription
prescription_trees = ttops_temp %>%
dplyr::ungroup() %>%
dplyr::left_join(
keep_trees
, by = "treeID"
) %>%
# tracking vars
dplyr::mutate(
# fill in keep tree flag
is_keep_tree = dplyr::coalesce(is_keep_tree, 0)
, orig_clump_n_trees_grp = clump_n_trees_grp
) %>%
dplyr::select(-c(clump_n_trees_grp))
# attach the new clumping to the trees
prescription_trees = prescription_trees %>%
dplyr::left_join(
# reclump
prescription_trees %>%
dplyr::filter(is_keep_tree == 1) %>%
st_clump_points(clump_dist_m = tree_clump_dist_m) %>%
sf::st_drop_geometry() %>%
dplyr::select(treeID, clump_n_trees_grp) %>%
dplyr::rename(new_clump_n_trees_grp = clump_n_trees_grp)
, by = "treeID"
) %>%
dplyr::mutate(
new_clump_n_trees_grp = forcats::fct_na_value_to_level(new_clump_n_trees_grp, level = "Cut tree")
)check the distribution of current vs prescription tree clump sizes
prescription_trees %>%
sf::st_drop_geometry() %>%
dplyr::count(
orig_clump_n_trees_grp, new_clump_n_trees_grp
) %>%
tidyr::pivot_wider(names_from = orig_clump_n_trees_grp, values_from = n, values_fill = 0) %>%
dplyr::arrange(new_clump_n_trees_grp) %>%
dplyr::rename(`New` = new_clump_n_trees_grp) %>%
kableExtra::kbl() %>%
kableExtra::kable_styling() %>%
kableExtra::add_header_above(
c(" "=1,"Old"=6)
)| New | Individual | 2-4 trees | 5-9 trees | 10-15 trees | 16-25 trees | >25 trees |
|---|---|---|---|---|---|---|
| Individual | 12 | 0 | 0 | 0 | 0 | 52 |
| 2-4 trees | 0 | 30 | 0 | 0 | 0 | 9 |
| 5-9 trees | 0 | 0 | 31 | 0 | 0 | 25 |
| 10-15 trees | 0 | 0 | 0 | 11 | 0 | 10 |
| 16-25 trees | 0 | 0 | 0 | 0 | 40 | 17 |
| Cut tree | 0 | 0 | 0 | 0 | 0 | 532 |
get the distance raster and openings vector data
ttops_temp2 = prescription_trees %>%
dplyr::filter(is_keep_tree == 1) %>%
dplyr::select(treeID, dbh_cm, tree_height_m, basal_area_m2) %>%
st_clump_points(clump_dist_m = tree_clump_dist_m)
clump_polys_temp2 = get_clump_summary(ttops_temp2)
dist_rast_temp2 = get_clump_dist_rast(clump_polys_temp2)
# what are the openings
dist_rast_temp2$openings_vect$openining_area_m2 %>% summary()## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.04 0.08 0.98 380.88 22.97 7460.28dist_rast_temp2$dist_rast %>% terra::aggregate(2, cores = 4) %>% terra::summary()## layer
## Min. : 0.0000
## 1st Qu.: 0.1707
## Median : 2.1610
## Mean : 3.2247
## 3rd Qu.: 4.9830
## Max. :23.4648plot the prescription
plt_fnl_temp2 =
ggplot() +
# distance
geom_tile(
data = dist_rast_temp2$dist_rast %>% terra::aggregate(2, cores = 4) %>% as.data.frame(xy = T) %>% rename(f=3)
, mapping = aes(x=x, y=y, fill = f)
) +
harrypotter::scale_fill_hp(
"mischief"
, na.value = "transparent"
, name = "distance to\nnearest tree (m)"
, limits = c(0,23)
) +
# openings
geom_sf(data = dist_rast_temp2$openings_vect, mapping = aes(color = openining_area_m2), fill = NA) +
scale_color_gradient(
low = "gray77", high = "gray11"
, labels = scales::comma_format(accuracy = 1)
, name = latex2exp::TeX("opening\narea ($\\m^2$)")
, limits = c(0,7500)
) +
# clumps
ggnewscale::new_scale_fill() +
geom_sf(
data = clump_polys_temp2
, mapping = aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
# tree points
geom_sf(data = ttops_temp2, color = "gray88", shape = ".") +
theme_void() +
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
# plot
plt_fnl_temp2
# save it
ggsave("../data/06_clumps_dist_presciption.jpg", width = 7, height = 7)highlight the openings
plt_open_temp2 =
ggplot() +
# clumps
geom_sf(
data = clump_polys_temp2
, mapping = aes(fill = clump_n_trees_grp)
, color = NA
) +
harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
# openings
ggnewscale::new_scale_fill() +
geom_sf(data = dist_rast_temp2$openings_vect, mapping = aes(fill = openining_area_m2), color = NA) +
scale_fill_gradient(
low = "gray77", high = "gray11"
, labels = scales::comma_format(accuracy = 1)
, name = latex2exp::TeX("opening\narea ($\\m^2$)")
, limits = c(0,7500)
) +
# tree points
geom_sf(data = ttops_temp2, color = "gray88", shape = ".") +
theme_void() +
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
plt_open_temp2
# save it
ggsave("../data/07_clumps_openings_presciption.jpg", width = 7, height = 7)combine plots with prescription
plt_fnl_temp2 + (plt_open_temp2 + theme(legend.position = "none")) &
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
ggsave("../data/08_clumps_combined_prescription.jpg", width = 10, height = 7)combine old vs new
(plt_fnl_temp + labs(subtitle = "Pre-Treatment")) +
(plt_fnl_temp2 + labs(subtitle = "Post-Treatment") + theme(legend.position = "none")) &
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
# save it
ggsave("../data/09_clumps_dist_compared.jpg", width = 10, height = 7)combine old vs new openings highlight
(plt_open_temp + labs(subtitle = "Pre-Treatment")) +
(plt_open_temp2 + labs(subtitle = "Post-Treatment") + theme(legend.position = "none")) &
theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))
ggsave("../data/10_clumps_openings_compared.jpg", width = 10, height = 7)table comparison
current vs target vs prescription
clump_n_trees_grp_summary_temp2 = get_clump_n_trees_grp_summary(
trees = ttops_temp2, clumps = clump_polys_temp2
)
clump_n_trees_grp_summary_temp2 %>% dplyr::glimpse()## Rows: 5
## Columns: 31
## $ stand_area_ha <dbl> 2, 2, 2, 2, 2
## $ clump_n_trees_grp <ord> Individual, 2-4 trees, 5-9 trees, 10-15 tr…
## $ n_trees <int> 64, 39, 56, 21, 57
## $ mean_dbh_cm <dbl> 48.15347, 44.26804, 45.02925, 38.01966, 43…
## $ mean_tree_height_m <dbl> 17.43778, 16.08613, 16.38036, 14.77738, 15…
## $ loreys_height_m <dbl> 22.07153, 23.74545, 22.43439, 21.57890, 23…
## $ basal_area_m2 <dbl> 14.270633, 8.535523, 11.883499, 3.273830, …
## $ clump_area_ha <dbl> 0.18087306, 0.08665170, 0.11212985, 0.0341…
## $ stand_n_clumps <int> 64, 16, 8, 2, 3
## $ trees_per_ha <dbl> 32.0, 19.5, 28.0, 10.5, 28.5
## $ basal_area_m2_per_ha <dbl> 7.135316, 4.267761, 5.941749, 1.636915, 5.…
## $ qmd_cm <dbl> 53.28276, 52.78832, 51.97963, 44.55264, 51…
## $ stand_trees_per_ha <dbl> 118.5, 118.5, 118.5, 118.5, 118.5
## $ stand_basal_area_m2 <dbl> 49.95235, 49.95235, 49.95235, 49.95235, 49…
## $ stand_basal_area_m2_per_ha <dbl> 24.97618, 24.97618, 24.97618, 24.97618, 24…
## $ pct_stand_basal_area <dbl> 0.28568490, 0.17087328, 0.23789668, 0.0655…
## $ pct_stand_n_trees <dbl> 0.27004219, 0.16455696, 0.23628692, 0.0886…
## $ stand_qmd_cm <dbl> 51.80347, 51.80347, 51.80347, 51.80347, 51…
## $ mean_dbh_in <dbl> 18.97247, 17.44161, 17.74152, 14.97975, 17…
## $ qmd_in <dbl> 20.99341, 20.79860, 20.47998, 17.55374, 20…
## $ stand_qmd_in <dbl> 20.41057, 20.41057, 20.41057, 20.41057, 20…
## $ mean_tree_height_ft <dbl> 57.19592, 52.76250, 53.72757, 48.46981, 52…
## $ loreys_height_ft <dbl> 72.39463, 77.88509, 73.58479, 70.77880, 75…
## $ basal_area_ft2_per_ac <dbl> 31.102844, 18.603171, 25.900086, 7.135313,…
## $ stand_basal_area_ft2_per_ac <dbl> 108.8712, 108.8712, 108.8712, 108.8712, 10…
## $ trees_per_ac <dbl> 12.9600, 7.8975, 11.3400, 4.2525, 11.5425
## $ stand_trees_per_ac <dbl> 47.9925, 47.9925, 47.9925, 47.9925, 47.9925
## $ stand_area_ac <dbl> 4.942, 4.942, 4.942, 4.942, 4.942
## $ clump_area_ac <dbl> 0.44693734, 0.21411636, 0.27707286, 0.0842…
## $ basal_area_ft2 <dbl> 153.60909, 91.87636, 127.91398, 35.23951, …
## $ stand_basal_area_ft2 <dbl> 537.6871, 537.6871, 537.6871, 537.6871, 53…target_data_temp %>%
dplyr::select(
clump_n_trees_grp
, tidyselect::starts_with("pct_stand_n_trees")
, tidyselect::starts_with("trees_per_acre_")
) %>%
dplyr::left_join(
clump_n_trees_grp_summary_temp2 %>%
dplyr::select(clump_n_trees_grp, pct_stand_n_trees, trees_per_ac) %>%
dplyr::rename(
pct_stand_n_trees_presc = pct_stand_n_trees
, trees_per_acre_presc = trees_per_ac
)
, by = "clump_n_trees_grp"
) %>%
dplyr::mutate(
dplyr::across(
tidyselect::where(is.numeric)
, ~ dplyr::coalesce(.x,0)
)
) %>%
dplyr::relocate(
pct_stand_n_trees_presc
, .after = pct_stand_n_trees_target
) %>%
dplyr::relocate(
trees_per_acre_presc
, .after = trees_per_acre_target
) %>%
dplyr::mutate(
dplyr::across(
tidyselect::starts_with("pct_stand_n_trees")
, ~ scales::percent(.x,accuracy = 1)
)
) %>%
kableExtra::kable(
caption = paste0(
"Current stand BA (ft2/ac): "
, clump_n_trees_grp_summary_temp$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand QMD (in): "
, clump_n_trees_grp_summary_temp$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
, "<br>Current stand TPA: "
, clump_n_trees_grp_summary_temp$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
)
, escape = F
, digits = 1
, col.names = c(
"", rep(c("current","target","prescription"),2)
)
) %>%
kableExtra::kable_styling() %>%
kableExtra::add_header_above(
c(" "=1,"% Trees"=3, "TPA"=3)
)| current | target | prescription | current | target | prescription | |
|---|---|---|---|---|---|---|
| Individual | 2% | 30% | 27% | 2.4 | 14.4 | 13.0 |
| 2-4 trees | 4% | 25% | 16% | 6.1 | 12.0 | 7.9 |
| 5-9 trees | 4% | 20% | 24% | 6.3 | 9.6 | 11.3 |
| 10-15 trees | 1% | 10% | 9% | 2.2 | 4.8 | 4.3 |
| 16-25 trees | 5% | 15% | 24% | 8.1 | 7.2 | 11.5 |
| >25 trees | 84% | 0% | 0% | 130.6 | 0.0 | 0.0 |
check the pre- and post-treatment summary
# call it
clump_n_trees_grp_summary_temp2 = get_clump_n_trees_grp_summary(
trees = ttops_temp2
, clumps = clump_polys_temp2
)
# table it
clump_n_trees_grp_summary_temp2 %>%
dplyr::mutate(presc = "post-treatment") %>%
dplyr::bind_rows(
clump_n_trees_grp_summary_temp %>%
dplyr::mutate(presc = "pre-treatment")
) %>%
dplyr::arrange(clump_n_trees_grp, desc(presc)) %>%
dplyr::select(
clump_n_trees_grp, presc
, n_trees
, mean_dbh_in
, qmd_in
, mean_tree_height_ft
, loreys_height_ft
, trees_per_ac
, basal_area_ft2_per_ac, pct_stand_basal_area, pct_stand_n_trees
) %>%
dplyr::mutate(
dplyr::across(
.cols = c(pct_stand_basal_area, pct_stand_n_trees)
, .fns = ~ scales::percent(.x, accuracy = 1)
)
) %>%
kableExtra::kbl(
digits = 1
, escape = F
, caption = "Pre- and Post-treatment overstory tree clump summary"
, col.names = c(
".", ""
, "trees"
, "mean<br>DBH (in)"
, "QMD (in)"
, "mean<br>Ht. (ft)"
, "Loreys<br>Ht. (ft)"
, "TPA"
, "BA<br>ft<sup>2</sup> ac<sup>-1</sup>"
, "%<br>stand BA"
, "%<br>stand trees"
)
) %>%
kableExtra::kable_styling() %>%
kableExtra::collapse_rows(columns = 1, valign = "top")| . | trees |
mean DBH (in) |
QMD (in) |
mean Ht. (ft) |
Loreys Ht. (ft) |
TPA |
BA ft2 ac-1 |
% stand BA |
% stand trees |
|
|---|---|---|---|---|---|---|---|---|---|---|
| Individual | pre-treatment | 12 | 24.3 | 26.6 | 66.2 | 80.2 | 2.4 | 9.3 | 4% | 2% |
| post-treatment | 64 | 19.0 | 21.0 | 57.2 | 72.4 | 13.0 | 31.1 | 29% | 27% | |
| 2-4 trees | pre-treatment | 30 | 17.8 | 21.8 | 52.8 | 81.4 | 6.1 | 15.8 | 7% | 4% |
| post-treatment | 39 | 17.4 | 20.8 | 52.8 | 77.9 | 7.9 | 18.6 | 17% | 16% | |
| 5-9 trees | pre-treatment | 31 | 13.5 | 15.7 | 45.6 | 65.7 | 6.3 | 8.4 | 4% | 4% |
| post-treatment | 56 | 17.7 | 20.5 | 53.7 | 73.6 | 11.3 | 25.9 | 24% | 24% | |
| 10-15 trees | pre-treatment | 11 | 11.5 | 12.4 | 41.9 | 50.7 | 2.2 | 1.9 | 1% | 1% |
| post-treatment | 21 | 15.0 | 17.6 | 48.5 | 70.8 | 4.3 | 7.1 | 7% | 9% | |
| 16-25 trees | pre-treatment | 40 | 18.3 | 21.4 | 54.6 | 76.7 | 8.1 | 20.2 | 9% | 5% |
| post-treatment | 57 | 17.2 | 20.4 | 52.4 | 75.9 | 11.5 | 26.1 | 24% | 24% | |
| >25 trees | pre-treatment | 645 | 13.6 | 15.3 | 46.2 | 61.6 | 130.6 | 166.3 | 75% | 84% |